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ABSTRACT 

Tsiganis et al. [Tsiganis, K., et al., 2005. Nature 435, 459-461] have proposed that 
the current orbital architecture of the outer solar system could have been established 
if it was initially compact and Jupiter and Saturn crossed the 2:1 orbital resonance 
by divergent migration. The crossing led to close encounters among the giant planets, 
but the orbital eccentricities and inclinations were damped to their current values by 
interactions with planetesimals. Brunini [Brunini, A., 2006. Nature 440, 1163-1165] 
has presented widely publicized numerical results showing that the close encounters led 
to the current obliquities of the giant planets. We present a simple analytic argument 
which shows that the change in the spin direction of a planet relative to an inertial 
frame during an encounter between the planets is very small and that the change in 
the obliquity (which is measured from the orbit normal) is due to the change in the 
orbital inclination. Since the inclinations are damped by planetesimal interactions on 
timescales much shorter than the timescales on which the spins precess due to the 
torques from the Sun, especially for Uranus and Neptune, the obliquities should return 
to small values if they are small before the encounters. We have performed simulations 
using the symplectic integrator SyMBA, modified to include spin evolution due to the 
torques from the Sun and mutual planetary interactions. Our numerical results are 
consistent with the analytic argument for no significant remnant obliquities. 



1. INTRODUCTION 



One of the fundamental questions of solar system formation is the origin of planetary spins. The 
spin axes of three of the four giant planets in our solar system are tilted significantly, with the present 
values of the obliquities (the angle between the spin axis and the orb it norma l ) of Ju piter, Saturn, 
Uranus, and Neptune equal to 3?1, 27°, 98°, and 30°, respectively. iBruninil (|2006al ) has recently 
proposed a unified mechanism for the origin of the ob liquities of the giant planets. The starting 
point of Brunini's model is the so-called Nice model ([Tsiganis et al.ll2005l ) for the establishment 
of the orbital architecture of the giant planets. In the Nice model, the outer solar system was 
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initially compact, with Jupiter and Saturn closer than their mutual 2:1 mean-motion resonance 
and all of the giant planets within ~ 18 AU of the Sun. The scattering of planetesimals initially in 
a disk starting just beyond the orbits of the planets caused Jupiter to migrate inward and Saturn, 
Uranus, and Neptune outward. When Jupiter and Saturn crossed the 2:1 mean-motion resonance, 
their orbital eccentricities were excited and the orbits of Uranus and Neptune were destabilized. A 
phase of close encounters among the giant planets followed, and Uranus and Neptune were scattered 
outward. The orbital eccentricities and inclinations of the giant planets were eve ntually damped 
to their present values by the interactions with planetesimals. iBruninil (|2006al ) has performed 
numerical simulations of the Nice model that follow the evolution of the planetary spin axes. He 
found that large obliquities are generated during the encounter phase and that the final obliquities 
are similar to the present values, including the 98° obliquity of Uranus. 



Bruninil (|2006bl ) has subsequently retracted his paper on the basis that simulations performed 



with another numerical code do not reproduce the previous results and show changes in the obliq- 
uities of only a few degrees in most runs. However, he did not provide any physical explanation as 
to why the later numerical results are more likely to be correct. Given the wide publicity of the 
first results, we thought it important to show in this paper a simple analytic argument (Section 2) 
that explains why significant obliquities of the outer planets cannot result from close encounters be- 
tween these planets in the Nice model. The analytic argument is veri fied by numerica l simulations 
(Section 3). We note t hat a s ummary of our results was presented bv lLee et al.l (|2006l ) prior to the 
retraction by IBruninil (|2006bl ). While our focus is on the Nice model, our analysis of the obliquity 
changes in encounters can also be adapte d to any scenario invol ving; close encounters among the 
planets, such as the scenario proposed by iGoldreich et al.l (|2004l . see also IChiang et al.l 12000 ) . in 
which several ice giants were formed in the outer solar system and all but two were ejected by 
encounters. 



2. ANALYTIC ARGUMENT 

Let us consider the change in the spin direction of a planet of mass M, radius R, and spin 
frequency uj, due to an encounter with a perturbing planet of mass M peTt . If we assume principle 
axis rotation, the spin angular momentum of the planet is 

L = (A + £)MR 2 ujs, (1) 

and the torque from the perturber at a distance r in the direction r is 

N = 3^%^(J 2 + q')MR 2 (f ■ s)(f x a), (2) 



where s is the unit vector of the spin direction, A is the moment of inertia of the planet normalized 
to MR 2 , and J2 is the quadrupole coefficient of the planet. The torque from the planet's oblate 
figure on a satellite locks the satellite to the planet's equator plane if the nodal precession period of 
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the satellite is sh ort compared to the timescale on which the spin direction of the planet is changed 
( Goldreichlll965l ). In the above equations, I is the orbital angular momentum (normalized to MR 2 u>) 
of the satellites whose orbits are coupled to the planet's equator, and q' / J% is t he ratio of the torque 
from the perturber on the satellites to that directly exerted on the planet (see IWard and Hamilton 



2004! for explicit definitions of I and q'). From dL/dt = N, the rate of change of the spin direction 
is 
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where Mq is the mass of the Sun and 
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(4) 



is the precession constant of the planet (plus its satellites) at the current orbital semimajor axis, 
A, of the planet. 

Because the torque decreases rapidly with the distance r, we expect the total change in the 
spin direction in an encounter, | As|, to be roughly equal to the rate of change at closest approach, 
| ds jdt | r , times the duration of closest approach, r p /v p : 
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where r p and v p are, respectively, the relative separation and speed at closest approach and the 
proportionality constant is a function of the encounter geometry. A more accurate expression for 
As can be obtained by neglecting solar gravity and using the two-body approximation during the 
encounter. Then 

+ axcco S (-l/e) ds 

-17 I "IT df, (6) 

- arccos(— 1/e) al \ " r / 

where e(> 1) and / are, respectively, the eccentricity and true anomaly of the relative orbit. In the 
coordinate system where the relative orbit is in the x-y plane and the closest approach is in the x 
direction, the integration in Eq. ([6]) yields 
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where 9 and <f> are the spherical coordinate angles specifying the orientation of the spin axis in 
the xyz system: s = sin # cos 05: + sin#sin</>y + cos#£. Contour plots of the magnitude, |As|, 
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as a function of 9 and 4> show that |As| reaches equal maximum values at four points with 9 = 
45° or 135° and <j) = 0° or 180° for any e > 1. For e = 1 (i.e., parabolic encounter), |As| = 
(ir/2)a(M per t/Mo)(A/r p ) 2 (A/vp)\ sin2#|, which is independent of 4> and maximum at 9 = 45° or 
135°. The plot of the maximum |As| for a given e as a function of e shows that it decreases 
monotonically with increasing e. Thus |Asj, as a function of e, 9, and 4>, is maximum for e = 1 
and 9 = 45° or 135°: 
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Since Uranus has the largest current obliquity and the largest changes in its spin direction 
come from encounters with Saturn, we consider first the change in the spin direction of Uranus due 
to an encounter with Saturn. As we shall see in Figs. [I] and [2 Saturn's orbit (of semimajor axis 
as) remains nearly circular, and an encounter between Saturn and Uranus must have Uranus on 
an eccentric orbit (of eccentricity eu), with Uranus usually near its perihelion. Thus the relative 
encounter speed at large separation is v w (GM /as) 1/,2 [(l + eu) 1 / 2 — 1] < 0.22(GM /os) 1 / 2 if 
e u ^ 0.5. Since the encounter phase occurs long after the formation of the regular satellites in the 
Nice model, the regular satellite systems must survive the encounter, and a reasonable lower limit 
to the separation r p at closest approach is the semimajor axis aTitan of the orbit of Titan around 
Saturn. For r p = a-pitam 2G(M S + M v )/r p > vl (where M s and M v are the masses of Saturn 
and Uranus, respectively), and v p = [v$ + 2G(M S + M u )/r p ] 1 / 2 w [2G(M S + M u )/r p ] 1 / 2 . So the 
maximum change in the spin direction of Uranus from an encounter with Saturn is 



l^ s lmax,U 



0?33(a Tit an/r p ) 3/2 , 



(9) 



if we adopt a\j = (4.6 Myr) 1 from iTremaind (|199ll ) 



For the same encounter between Uranus and Saturn, the maximum change in the spin di- 
rection of Saturn is |As|U axS = 0?ll fa T ^ar 1 /r> ) 3/2 if as = (0.26 Myr)" 1 (^remaind fl99lh . If we 
also adopt on = (17.6 Myrj" 1 from Tremaing ( 199 ll ). the changes from an encounter between 
Neptune and Saturn are nearly identical to those from an encounter between Uranus and Saturn, 
with |As|^ ax N w |As|^ ax U and lAs^g ^ lAa^g, because Uranus and Neptune have sim- 
ilar masses and au^4fj ~ a N^.j\j (see Eq. (|8|)). For an encounter between Uranus and Neptune, 
|As|^ axU m |As|^ axN 0?3(aoberon/ r p) 3//2 - We do not consider encounters with Jupiter, since 
such encounters usually result i n systems that do no t resemble the solar system by, e.g., ejecting 
the planet encountering Jupiter ( Tsiganis et al.1 2005 ). 



The precession constants adopted above include the contributions from the satellites of the 
planets. In the case of Uranus, the contributions are mainly from Titania and Oberon. However, 
the nodal precession periods of Titania and Oberon (> 1300 yr due to the obla teness of Uranus alone 
and ~ 200 yr with the inclusion of the secular interactions among the satellites; lLaskar and Jacobson 
19871 ) are much longer than the encounter timescale with either Saturn (r p /v p ~ 2 days for r p ~ 
«Titan) or Neptune (r p /v p ~ lday for r p ~ a ohprnn), and thei r orbits would not be able to follow 
the equator of Uranus during the encounter (IGoldreichl 1 19651 ) . This would reduce a\j and hence 
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the already small maximum change in the spin direction of Uranus by a factor of 5.3. Similarly, 
without the contributions from the satellites to the precession constants, the maximum changes 
in the spin directions of Neptune and Saturn are reduced by a factor of 5.7 and 3.7, respectively. 
Then the maximum change in the spin direction is < 1° for any pair-wise encounter among Saturn, 
Uranus, and Neptune, even if we consider an extremely close encounter with r p equal to twice the 
sum of the planetary radii, which the regular satellite systems are unlikely to survive. 

We have just shown that the change in the spin direction in inertial space during an encounter 
between the planets is very small. Thus any change in the obliquity (which is defined as the angle 
between the spin axis and the orbit normal) must be due to a change in the orbital inclination. As 
we shall see in, e.g., Figs. Q] and [21 the orbital inclinations of the ice giants can reach maximum 
values of ~ 10°-15° during the encounter phase. However, because the timescale on which the 
inclination is changed by encounters and damped by planetesimals (< 1 Myr) is much less than the 
precession period of the spin axis about the orbit normal due to solar torque (> 27r/au = 29 Myr 
for Uranus), there is not enough time for the spin axis to precess much due to solar torque, and the 
spin axis does not stray far from its initial position nearly perpendicular to the initial orbital planes 
of the planets. Therefore, we expect the obliquities to return to small values as the inclinations are 
damped to their present values. 



NUMERICAL SIMULATIONS 



To verify the analytic argument in Section 2 and to eliminate the possibility of unexpected 
secular spin-orbit effects, we have performed numeri cal simulations of t he Sun and the four giant 
planets using the symplectic iV-body code SyMBA (IDuncan et al.lll998l ). which we have modified 
to include spin evolution. For the orbital evolution, we ignore the negligible effects of the oblateness 
of the planets, and we use impo sed migration and da mping of eccentricity and inclination to model 
the effects of the planetesimals. iLee and Peald (|2002l ) have described how SyMBA can be modified 
to include forced migration and damping of eccentricity and inclination. The unit vector s k of the 
spin direction of planet k is evolved according to 



ds k 
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(10) 



where the symbols are as in Eq. ([3]), but with the subscript for the Sun and j,k ^ for the 
planets. In Eq. (|10p . the first term is due to the torque from the Sun and the second term is due 
to the torques from the other planets. The recursively subdivided time step used by SyMBA to 
handle close encounters between the planets is also implemented for the spin evolution due to the 
planetary torques. In the Appendix, we describe in more details how SyMBA is modified to include 
spin evolution and the tests that were performed. 
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We present the results from 4 series of runs, each with 30 simulations. The simulations were 
performed with a time step of 0.125 yr. For the precession constants at the current orbital semimajor 
axes, we adopt otj = (0.077 Myr) -1 for Jupiter , as = (0.26 Myr) -1 for Saturn, = a\j = 
(4.6 Myr) -1 for the initially inner ice giant, and a\ 2 = on = (17.6 M yr) -1 for the in itially outer 
ice giant. These values include the contributions from the satellites ( Tr emainel 1 1 99 ll ) . Although 



the inner and outer ice giants can switch positions (see, e.g., Fig. [I]), the fact that we adopt the 
parameters of Uranus (Neptune) for the initially inner (outer) ice giant is not important, because 
Uranus and Neptune have similar masses and aijAfj ~ un^. The four series of runs have the 
following initial conditions in common: the orbital eccentricities are 10 -3 ; both the orbit normals 
and the spin axes are tilted by 10 -3 radians (or 0?057) from the z-axis of the inertial frame; the 
mean anomalies, the arguments of perihelion, and the longitudes of the ascending nodes on the x-y 
plane of both the orbits and the equators are random variables. 

In series I, the initial orbital semimajor axes and t he imposed mig ration and damping of 



eccentricity and inclination are similar to those adopted by iBruninil (|2006al ). The initial semimajor 
axis of Jupiter is aj = 5.45 AU. The initial semimajor axis of Saturn is varied in the range 8.38- 
8.48 AU, whereas those of the ice giants are varied in the ranges 9.9-12 AU and 13.4-17.1 AU, while 
keeping the initial orbital separation of the ice giants larger than 2AU. The imposed migration 
rate is ^ 

a k = exp(-i/r), (11) 

T 

where Aa k is the difference between the current semimajor axis, A k , and the initial semimajor axis. 
The migration timescale r is varied between 1 and 10 Myr, and the total time span of a simulation 
is lOr. The imposed eccentricity and inclination damping rates are 

ejt = -e fc /(2r e ), i k = -i k /{2T i ) 1 (12) 

respectively, with r e = Tj = r/10. The eccentricity and inclination dampings are applied only if 
the orbital distance of a planet is between just outside the initial orbit of the outer ice giant and 
30 AU. Series II is similar to series I, but with r e = Tj = r/20. 



Series III and IV have initial orbital semimajor axes closer to those adopted by lTsiganis et al 



(|2005l ) and use a different prescription for migration and damping. The initial aj = 5.45 AU, and 
the initial semimajor axes of Saturn and the inner and outer ice giants are varied in the ranges 
8.0-8.5 AU, 11-13 AU, and 13.5-17 AU, respectively, with the ice giants at least 2AU apart. The 
main difference from series I and II is that the semimajor axis of the inner ice giant is about 1 AU 
larger. For the imposed migration, we use Eq. (jlip for Jupiter and Saturn and 

a k = ^ (~) expH 2 /(2r 2 )] (13) 



for the ice giants. With the inner ice giant further from Saturn initially, the migration rate in 
Eq. (fl~3j) ensures that the ice giants do not migrate too far by the time of the Jupiter-Saturn 2:1 
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resonance crossing, as seen in the simulations of lTsiganis et al.l (|2005l ). For all planets, the damping 
rates are 

e-k/ek = -K e \a k /a k \, i k /i k = -Ki\a k /a k \, (14) 
and they are app lied at all times. T his form of damping has been used to model planet-disk 



interactions (e.g., iLee and Peald 120021 ) . and it has the reasonable property of tying the damping 
rates to the migration rate (i.e., less damping when the migration slows down). We adopt K e = 
Ki = 10 for series III and K e = Ki = 5 for series IV. 

The fractions of cases with all four giant planets surviving on reasonable orbits are 47%, 73%, 
57 %, and 50% in series I-IV, respectively. They are all roughly consistent with the 67% reported 



bv lTsiganis et al.l (|2005l ). given the small number of simulations in Tsiganis et al. (43 simulations) 
and in each of our series (30 simulations). The results from two typical simulations with all four 
giant planets surviving on reasonable orbits — one from series I and the other from series III 
- are shown in Figs. Q] and [21 Fig. [T] is an example where the ice giants switch positions. In 
each Figure, panel a shows the evolution of the orbital semimajor axes a, the perihelion distances 
q = a(l — e), and the aphelion distances Q = a(l + e), while panels b, c, and d show the evolution 
of the orbital inclinations i, the tilts of the spin axes from the z-axis of the inertial frame (which is 
nearly perpendicular to the initial orbital planes of the planets), and the obliquities e, respectively. 
It can be seen from the correlation of the changes in panels b and d of each Figure that large 
changes in the obliquities of Saturn and the ice giants during close encounters are due to changes in 
orbital inclinations. The changes in the directions of the spin axes relative to inertial space during 
the encounters are minuscule (panel c). As the inclinations are damped to small values, the spin 
axes keep their orientations in inertial space, and the obliquities are similarly decreased. These 
behaviors are exactly what we expect from the analytic argument in Section 2. 

For Jupiter, we note that both the obliquity and the tilt of the spin axis from the inertial 
z-axis increase to 1-3 degrees in Figs. Q] and [2j This could be due to secular spin-orbit coupling, 
because the spin precession due to solar torque is quite fast for Jupiter (aj = (0.077Myr) _1 ), but 
a detailed examination of this process is beyond the scope of the present paper. 

Among the cases with all four giant planets surviving on reasonable orbits (68 out of 120 
simulations in all four series), there are only four cases where one or both of the ice giants have 
final obliquities greater than 5°, including one case where one of the ice giants (Neptune) has a 
final obliquity greater than 10° (en = 27°). However, these cases involve one or two extremely close 
encounters between the planets (usually between Saturn and the initially inner ice giant), for which 
the regular satellite systems are unlikely to survive. Even if the regular satellites could survive 
these encounters, their orbits would not be able to follow the equators of the planets during the 
encounters, and the final obliquities obtained in our simulations (which include the contributions 
from the satellites) should be reduced by a factor of 5-6 (see Section 2) and < 5° even in the most 
extreme case. 
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CONCLUSIONS 



We have analyz ed the evolution of the obliquities of the giant planets in the Nice model 
(jTsiganis et al.l 120051 ). where the Jupiter-Saturn 2:1 resonance crossing led to close encounters 
among the giant planets and the orbital eccentricities and inclinations were eventually damped 
to their present values by the interactions with planetesimals. We presented a simple analytic 
argument which shows that the change in the spin direction of a planet in inertial space during an 
encounter between the planets is very small and that the change in the obliquity (which is measured 
from the orbit normal) is due to the change in the orbital inclination. Since the inclinations 
are changed by encounters and damped by planetesimal interactions on timescales much shorter 
than the timescales on which the spins precess due to solar torques, especially for Uranus and 
Neptune, the obliquities return to small values if they are small before the encounters. We have 
performed simulations using the symplectic integrator SyMBA modified to include spin evolution. 
The numerical results are consistent with the analytic argument leading to no significant remnant 
obliquities. The analytic arg ument provides a physical explanation as to why the numerical results 
reported bv lBruninil (2006a). which are not confirmed by later simulations by him using another 
numerical code (jBruninil l2006bl ) or the simulations in this paper, are indeed incorrect. 



Giant impacts in the late stages of the formation of Uranus and Neptune remain a plausible 
explanation for the obliquities of the ice giants, but a giant impact o rigin for Saturn's obliquity is 



problematic because of its large mass and spin angular momentum (jLissauer and Safronovlll991 



Pones and Tremaindll993l ). Alternatively, the obliquities of the giant planets could be generated 
by any proces s that twists th e total angular momentum vector of the solar system on a timescale 
of ~ 0.5 Myr (|Tr emaind 1 1 99 ll ) . However, the similarity between Saturn's spin-axis precession rate 



and the nodal regression rate of Neptune's orbit, and the near match of Saturn's 27° obliquity to 
that of the Cassini state 2 of the secular spin-orbit resonance, strongly argue for the capture of 
Saturn into the Cassini sta t e whose obliquity increases during the dispersal of the planetesimal disk 



(Ward and Hamilton 



2004 



Hamilton and Ward! 12004 ) . 
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A. NUMERICAL METHODS 



In this appendix we describe how the symplectic integrator SyMBA (jDuncan et al.lll998l ) is 
modified to include the evolution of the spin axes according to Eq. (I10p and the tests that were 
performed. Our algorithm differs from the Lie-Poisson integrator for rigid body dynamics by 



Touma and Wisdom! (|1994l ) in assuming principle axis rotation and handling close encounters. 
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SyMBA is based on a variant of the I Wisdom and Holmanl (|199ll ) method, with the gravitational 



A^-body Hamiltonian written in terms of positions relative to the Sun and barycentric momenta, 
and employs a multiple time step technique to handle close encounters. The potential of the 
gravitational interaction between planets j and k is of the form —Gm^m^/rj^, which we simplify 
as —1/r in the following discussion. In SyMBA, a set of cutoff radii Ri > Ri > ■ ■ ■ is chosen, 
with R\ a few times the mutual Hill radius, and the potential —1/r is decomposed into YleLo Vli 
or equivalently the force F into Y^=o^t = ~Y^Lo^^l^ r = ~ Y^tLo Pe{r)r/r 3 . The properties 
of the partition functions Pg(r) are such that (i) Pi (except Po) increases smoothly from zero at 
r < Ri + 2 to one at r = Ri + \ and then decreases smoothly to zero at r > Ri and (ii) Po increases 
smoothly from zero at r < R<i to one at r > R\. The force F g_ is to be applied with a time step hi, 
with hi/he+i an integer (= 3 in the standard implementation). A single step of overall time step 
ho of the SyMBA algorithm can be represented by the following sequence of substeps: 

-E'sun ( -jr J -Elnt.o ( j E recur (h )E int fi ( ) Es un ( ) . (Al) 



In the substep £sun(^o/2), each planet takes a linear drift in its heliocentric position (corresponding 
to a linear drift in the position of the Sun from the barycenter) for time ho/ 2. In the substep 
Pint, o(^o/2), each planet receives a kick to its linear momentum due to the i = level forces from 
the other planets for time ho/2. In the recursive substep P r ecur(^o)i a planet that is not having 
close encounters (i.e., a planet whose separation from any other planet is greater than R\) evolves 
along a Kepler orbit for time ho, while a pair of planets having an encounter have their time steps 
recursively subdivided to a maximum level of ^ max — with the force at lev el I, Ff, applied wit h a 



time step hi — if their separation r > P^ max +i during the overall step (see lDuncan et al.lll998l for 
details). 

In Eq. (jlOp . the torque on the spin axis of planet k from planet j is of the form 2a/ c (M ? /Mo) x 
(Ak/rjkY{r jk ■ s^.)(f jfc x s^), which we simplify as T = r _3 (f • s)(r x s) in this discussion. As 
for the force F, we decompose T into Y^iLo^e = -^H r ) r ~ 3 (^ • s)(f x s) using the partition 

functions Pi{r), and Ti is to be applied with a time step h(_. Thus our modified algorithm (with 
imposed migration and damping) can be represented by the following sequence of substeps: 

Ea (~2^) (if) ESUn ^un f "j) E ' intfi (if) ^ecur(^o) 

*k„ (f ) « (f) *- (f ) (f ) * (f ) ■ (A2) 

In the substep E a {ho/2), the orbital semimajor axis a of each planet is changed according to the 
exact sol ution to the imposed migration rate (with all of the other orbital elements fixed) for time 



ho/ 2 (see iLee and Peald 12002 for details). Similarly, the orbital eccentricity e and inclination i of 
each planet are changed according to the exact solutions to the imposed damping rates for time 
ho/2 in the substep E e i(ho/2). In the substep E^^(ho/2), the spin of each planet is evolved due 
to the solar torque (i.e., the first term in Eq. (|10p ) alone for time ho/2, which is a rotation of 
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Sfc by an angle 2afc(^4fc/ro/c) 3 (^ofc • s k)ho/2 about r fc, because r fc is fixed during this substep. 
The substep E' int (ho/2) is the same as £ , ; n t i o(^o/2) in the original SyMBA algorithm, but it also 
evolves the spin direction of each planet due to the I = level torques from the other planets for 
time ho/2. Unlike the substep applying the solar torque alone, there is not an exact solution for 
the change in the spin direction with the torques from multiple planets, and we use the midpoint 
method. Similarly, ^ecur(^o) is the same as E recnr (ho), but it also applies the torque at level £, Ti, 
with a time step hi, recursively down to the maximum level £ mSuX , for a pair of planets having an 
encounter. 

We have tested our implementation of the spin evolution with several tests. In the first set of 
tests, we checked that the code produces the correct spin axis precession around the orbit normal 
in the presence of the solar torque alone. We integrated Jupiter alone with a variety of orbital 
eccentricity e (up to 0.2). We also integrated the four giant planets with their current e and i, but 
with their masses reduced by a factor of 10~ 10 (so that they are on unperturbed Kepler orbits). 
A variety of obliquities e were used, and we integrated for 10 yr using a time step of 0.25 yr. We 
found a small variation in e (up to ~ 10~ 3 degrees) over the period of an orbit, which is expected 
since Eq. (fT0|) is not orbit-averaged, but there was no detectable secular change in e. The average 
spin precession rate (measured by the total change in the longitude of the ascending node of the 
equator on the orbital plane) agreed with the analytic result a(l — e 2 ) -3 / 2 cose to better than one 
part in 5 x 10 . 

In the second set of tests, we checked that the code produces the correct capture into (and 
es cape from) secular spin-orb it resonance. The simulations were simplified versions of the ones 
in lHamilton and Ward! (|2004l ). We integrated Saturn alone with a time step of 0.58 yr. Initially, 
o = 9.54 AU, e = 10 -3 and i = 0?064 (the last being the amplitude associated with the per- 
turbation to Saturn's orbital plane due to the nodal regression of Neptune). We imposed an 
orbital nodal precession rate of g = — a(0.89 + 0.23e - */ T ) (magnitude decreasing with time) or 
g = — a(1.12 — 0.23e~'/ r ) (magnitude increasing with time), with r = 2.5-10 x 10 8 yr. In the cases 
where \g\ was decreased from above a to below a, Satur n's spin was captured into C assini state 2, 
with the obliquity librating about the analytic result in IWard and Hamiltonl (J200J). In the cases 
where \g\ was increased from below a to above a, the sp in was first captured into Cassini state 
1, with the obliquity librating about the analytic result in lWard and Hamiltonl (120041 ) . There was 
a jump in the obliquity when Cassini state 1 disappeared, and the obliquity subsequently oscil- 
lated about a constant valu e. T he results agreed with the a nalytic theory and numerical results of 
Ward and Hamiltonl (|200J) and lHamilton and Ward! (|200J). 



As a very stringent test of the spin evolution due to close mutu al planetary in t eract ions, we 
performed simulations similar to the binary planet test described in iDuncan et al.l (119981 ) . where 
two planets are in a bound orbit and their center of mass orbits a star. The center of mass of the 
binary was placed in a 5.2 AU circular orbit about a solar-mass star. The initial semimajor axis 
of the relative orbit of the binary was 0.0125 AU and the initial eccentricity was 0.4. One planet 
had a mass of M 1 = 0.9 x 1O~ 3 M and the other M 2 = 1.1 x 1O~ 3 M . Both planets had the 
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same precession constant as Jupiter, and their initial obliquities were E\ = 60° and £2 = 30°. We 
integrated this system for 500 yr (or about 16000 orbital periods of the binary). Because the torques 
of the planets on each other dominated over the torques from the Sun, the average precession rates 
of the spin axes around the normal to the relative orbit are known analytically, if we ignore the 
solar torques. We found that the numerical precession rates agreed with the analytic results to 
within 0.0008 and 0.0022 fractionally for planets 1 and 2, respectively. The small discrepancies 
were not sensitive to the time step (we used 0.25 yr and smaller), and they are likely real and due 
to solar perturbations. Solar perturbations also caused the eccentricity of the relative orbit to vary 
between 0.4 and 0.3964 and the obliquities to show small variations on the same timescale. 
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Fig. 1. — Evolution of the orbits and spin axes of the giant planets in a simulation from series I 
(see text for the details of the setup), (a) Orbital semimajor axis (a) and the minimum (g) and 
maximum (Q) heliocentric distances. Note that the ice giants switch positions in this simulation, 
(b) Orbital inclination, (c) Tilt of the spin axis in inertial space, (d) Obliquity. The planets are 
Jupiter (red), Saturn (green), Uranus (blue), and Neptune (magenta). The orbital inclination and 
the tilt of the spin axis in inertial space are measured relative to the z-axis of the inertial frame, 
which is nearly perpendicular to the initial orbital planes of the planets, while the obliquity is the 
angle between the spin axis and the orbit normal. The directions of the spin axes show very little 
change in inertial space. The obliquities reflect primarily the inclinations and return to small values 
as the inclinations are damped. 




Fig. 2. — Same as Fig. [TJ but for a simulation from series III. The ice giants do not switch positions 
in this simulation. 



